#include <stdio.h>

#include "Jphi.h"
#include "Constants.h"
#include <math.h>
#include "Parameters.h"
#include "Grid.h"
#include "DynamicArray.h"

using namespace std;

extern Parameters* ppmts;


Jphi::Jphi(Grid* g) {
    grid = g;
    num_r = grid->getNumberR();
    num_z = grid->getNumberZ();
    left_r = grid->getLeftR();
    down_z = grid->getDownZ();
    h_r = grid->getHR();
    h_z = grid->getHZ();
    delta_s = grid->getDelta();
    alpha = ppmts->getAlpha();
    beta_pl = ppmts->getBetaPl();
    psi_g =ppmts->getPsiG();
    r0 = ppmts->getR0();
    I  = ppmts->getI();
    pl_left_r = ppmts->getPlasmaLeftR();
    pl_right_r = ppmts->getPlasmaRightR();
    pl_up_z = ppmts->getPlasmaUpZ();
    pl_down_z = ppmts->getPlasmaDownZ();

    jphi_0 = AllocateDynamicArray<float>(num_r, num_z);
    jphi_m = AllocateDynamicArray<float>(num_r, num_z);
    jphi_tmp = AllocateDynamicArray<float>(num_r, num_z);
    jphi_s = AllocateDynamicArray<float>(num_r, num_z);
};

Jphi::~Jphi(){

    FreeDynamicArray<float>(jphi_0);
    FreeDynamicArray<float>(jphi_m);
    FreeDynamicArray<float>(jphi_tmp);

    FreeDynamicArray<float>(jphi_s);
};


float** Jphi::on_grid(float** psi_grid, float** jphi) {
    float r, z = 0;
    float sum = 0;
    float gamma = 0;
    //innormalized jphi calculation
    for (int j = 0; j < num_z; j++) {
        z = down_z + j*h_z;
        for (int i = 0; i < num_r; i++) {
            r = left_r + i*h_r;
                 //check for Omega area (ext current excluded)
           
            if (psi_grid[i][j] > psi_g)
            {
                jphi[i][j] = pow((psi_grid[i][j] - psi_g), alpha)*
                        ((r / r0) * beta_pl + (r0 / r)*
                        (1 - beta_pl));
                sum = sum + jphi[i][j];
            }
            else {
                jphi[i][j] = 0;
            }
        }
     }
    sum = sum * delta_s;
    gamma = I/sum;
    printf("gamma: %f", gamma);
    for (int j = 0; j < num_z; j++) {
        for (int i = 0; i < num_r; i++) {
            jphi[i][j] = jphi[i][j] * gamma;
        }
    }
    return jphi;

}

float** Jphi::on_grid_initial() {
    for (int j = 0; j < num_z; j++) {
        for (int i = 0; i < num_r; i++) {
            if ( (i == num_r/2) && (j == num_z/2) ) {
                jphi_0[i][j] = I/delta_s;
   
            }
            else {
                jphi_0[i][j] = 0;
            }

        }
    }
    return jphi_0;

};


float Jphi::get_point(float r, float psi_) {
    if (psi_ > psi_g) {    
    return pow((psi_ - psi_g), alpha)*
            ((r / r0) * beta_pl + (r0 / r)*
            (1 - beta_pl));}
    else return 0;
   
};


